自由放置的刚体在地震下的摆动响应

自由放置的刚体在地震下的摆动响应1. 运动方程2. 摆动周期3. 能量耗散4. 在余弦激励下的响应4. 在地震激励下的响应

1. 运动方程

$\ddot u_g(t)$ $o$ $o'$ 发生摆动。

该刚体的运动方程为：

I_O\ddot\theta+mgRsin(-\alpha-\theta)=-m\ddot u_gRcos(-\alpha-\theta),\theta<0\\ I_O\ddot\theta+mgRsin(\alpha-\theta)=-m\ddot u_gRcos(\alpha-\theta),\theta>0 \quad

$I_O=(4/3)mR^2$ $I_O$ ，得：

\ddot\theta(t)=-p^2\{sin(\alpha sgn[\theta(t)]-\theta(t))+\frac{\ddot u_g}{g}cos(\alpha sgn[\theta(t)]-\theta(t))\}

$p=\sqrt{(3g/4R)}$ $rad/s$ $sgn$ 为符号函数。

2. 摆动周期

$\theta_0$ $\theta_0$ 的的状态摆回到竖直状态过程中的运动方程为：

I_O\frac{d^2\theta}{dt^2}=-mgRsin(\alpha-\theta)

$\alpha<20^\circ$ ），上述方程可以简化为线性形式：

I_O\ddot\theta-mgR\theta=-mgR\alpha

$I_O$ ,得到：

\ddot\theta-p^2\theta=-p^2\alpha

$\theta=\theta_0,\dot\theta=0$ ，该方程的解为：

\theta=\alpha-(\alpha-\theta_0)cosh(pt)

$t=T/4$ $\theta=0$ ，带入上式：

0=\alpha-(\alpha-\theta_0)cosh(p\frac{T}{4})

即

cosh\frac{pT}{4}=\frac{1}{1-\theta_0/\alpha}

则其周期为：

T=\frac{4}{p}cosh^{-1}(\frac{1}{1-\theta_0/\alpha})

上述方程绘制出来如下图所示。


4
1
theta_alpha=0:0.02:0.98;
2
Tp_4=acosh(1./(1-theta_alpha));
3
plot(theta_alpha,Tp_4);
4
grid on;xlabel('\theta_0/\alpha');ylabel('Tp/4');

$\alpha$ 时，其摆动周期接近无限大。物理意义上来说，这种状态下，刚体的重心在地面的投影落在支撑点上，处于平衡状态，不发生转动。一旦有所扰动，要么往回摆动，要么倾翻。

3. 能量耗散

刚体在摇摆的方向发生变化时，会与地面发生一次冲击，根据角动量守恒，可以得到：

I_O\dot\theta_1-2mRb\dot\theta sin\alpha=I_O\dot\theta_2

$\dot\theta^2$ 来描述，由于与地面碰撞产生能量损失，假设碰撞后刚体的动能为碰撞前的动能乘以一个系数r，即：

\dot\theta_{n+1}^2=r\cdot\dot\theta_n^2

结合以上两式，可得：

r=[1-\frac{3}{2}sin^2(\alpha)]^2

4. 在余弦激励下的响应

$p=2,\alpha=15^\circ$ 。则不同脉冲幅值下，刚体的转角和角速度分别如下图。

$\alpha$ ，即刚体的重心在地面的投影已经落到了支点的外侧，发生了倾翻。

为了观察第二节描述的刚体摆动周期随其转角无限增大的情况，这里可以逐步增大脉冲幅值至临界值，如下图所示。随着脉冲幅值的微幅增大，刚体的角速度接近于0的时间越长，在这段时间里，刚体近乎于静止状态。


61
1
function nonlinearrock(alpha,p)
2
    p=2;alpha=deg2rad(15);
3
    ap=0.315034034;
4
    omega=pi;    tstop=2;
5
    tg=0:0.1:10; tcos=0:0.1:tstop;
6
    ag=zeros(1,length(tg));
7
    ag(1:length(tcos))=ap*cos(omega*tcos);
8
    subplot(3,3,3)
9
    plot(tg,ag); xlabel('time(s)');ylabel('脉冲幅值(g)');title(['A=',num2str(ap,9),'g']);grid on;
10
11
    %%%
12
    tspan=[0,10];
13
    inity=[0,0];
14
    option=odeset('RelTol',1e-8,'AbsTol',1e-10,'Events',@EventFcn);
15
    while tspan(1)<tspan(end)
16
        [rt_temp,ry_temp,te,ye,ie]=ode45(@wffc,tspan,inity,option);
17
        subplot(3,3,6);plot(rt_temp,ry_temp(:,1)/alpha,'black');hold on;
18
        subplot(3,3,9);plot(rt_temp,ry_temp(:,2),'r'); hold on;
19
        if isempty(ie)
20
            break;
21
        end
22
        if ie(end)==1
23
            inity=ye(end,:);
24
            inity(2)=inity(2)*(1-1.5*(sin(alpha))^2);
25
            tspan=[te(end),tspan(end)];
26
            if abs(inity(2))<1e-6
27
                break;
28
            end
29
            continue;
30
        end
31
        if ie(end)==2
32
            break;
33
        end
34
    end
35
    subplot(3,3,6);xlabel('time(s)');ylabel('\theta/\alpha');hold off;grid on;
36
    subplot(3,3,9);xlabel('time(s)');ylabel('\theta''');hold off;grid on;
37
    
38
    function dy=wffc(t,y)
39
        dy=zeros(2,1);
40
        sgn=sign(y(1));
41
        if sgn==0
42
            sgn=sign(y(2));
43
        end
44
        dy(1)=y(2);
45
        if t<=tstop
46
            dy(2)=-p^2*(sin(alpha*sgn-y(1))+ap*cos(omega*t)*cos(alpha*sgn-y(1)));
47
        else
48
            dy(2)=-p^2*sin(alpha*sgn-y(1));
49
        end
50
    end
51
52
    function [value,isterminal,direction]=EventFcn(~,y)
53
        value(1)=y(1);
54
        isterminal(1)=1;
55
        direction(1)=0;
56
        
57
        value(2)=abs(y(1))-alpha;
58
        isterminal(2)=1;
59
        direction(2)=1;
60
    end
61
end

4. 在地震激励下的响应


97
1
function [rt,ry]=rock1(alpha,T,tg,ag)
2
% 求解参数为alpha、T的刚体在时间和幅值依次为tg、ag的地震作用下的响应。
3
% tg\ag:地震波的时间列和幅值列,ag的单位为g。
4
p=2*pi/T;
5
rt=tg;                      % 存储时间点
6
ry=zeros(length(tg),2);     % 用于存储每个时刻的响应结果，含有两列，第一列为角度，第二列为角速度
7
% 静止状态下，只有当地震加速度的幅值的绝对值abs(ag)>tan(alpha)时，刚体才开始摆动
8
% 寻找地震波中第一个幅值的绝对值超过tan(alpha)的点的下标index。
9
index=find(abs(ag)>tan(alpha),1); 
10
if isempty(index)
11
    return;
12
end
13
sgn=sign(ag(index));    %该值是正还是负
14
if index==1             %如果第一个值就足够大，则开始计算的时刻则为地震波的第一个时刻 
15
    tgtrig=tg(1);
16
    index=2;
17
else                    %否则，通过插值求得地震幅值等于tan(alpha)的时刻tgtrig.
18
    tgtrig=interp1([ag(index-1),ag(index)],[tg(index-1),tg(index)],sgn*tan(alpha));
19
end
20
    
21
%% 求解
22
%option:ode45的求解参数设置，RelTol为相对误差，AbsTol为绝对误差。
23
%Events参数指向一个函数，当函数的值为0时，则触发事件。该函数有两个作用：
24
% 1.当刚体摇摆方向切换时，中断求解，并将摇摆的角速度乘以一个比例系数。
25
% 2.当刚体的摇摆角度超过alpha时，认为刚体侵翻，终止计算。
26
option=odeset('MaxStep',0.025,'RelTol',1e-8,'AbsTol',1e-10,'Events',@EventFcn);
27
inity=[0,0];                % 计算的初始值
28
tstart=tgtrig;              % 开始求解的时刻，该时刻起，刚体才开始摆动
29
while tstart<tg(end)        % 每一步循环求解，每次只求解一个时间段，直到地震波的最后一刻
30
    %调用ode45进行求解，求解的时间段为[tstart,tg(index)]，该时间段求解完后，index加1，然后循环求解
31
    % @wffc为待求解的微分方程，该方程有3种可能的终止结果
32
    % 1：正常结束，完成时间段[tstart,tg(index)]段的求解
33
    % 2: 由于摇摆方向改变导致终止
34
    % 3：由于摇摆角度超过alpha导致终止。
35
    % rt_temp，ry_temp为求解后的时间、幅值
36
    % te为触发事件时的时间，ye为触发事件时的响应，
37
    % ie的值指示终止的原因，如果ie=1，则是上方可能2，如果ie=2，则是上方可能3
38
    [rt_temp,ry_temp,te,ye,ie]=ode45(@wffc,[tstart,tg(index)],inity,option);
39
    subplot(3,1,2);plot(rt_temp,ry_temp(:,1)/alpha,'r');hold on; %
40
    subplot(3,1,3);plot(rt_temp,ry_temp(:,2),'r');hold on;
41
    if rt_temp(end)==tg(index)      % ode45正常终止
42
        ry(index,:)=ry_temp(end,:); % 保存tg(index)时刻的响应
43
        inity=ry_temp(end,:);       % 该时间段最后时刻的结果作为下一时间段求解的初始值
44
        tstart=tg(index);           % 下一时间段的起始时刻为tg(index)
45
        index=index+1;              % 循环开始下一时间段的求解
46
        continue;
47
    end
48
    if ie(end)==1                   % 如果由于摇摆方向改变导致终止
49
        inity=ye(end,:);            % 下一次求解的初始值为中断求解时的响应
50
        inity(2)=inity(2)*(1-1.5*(sin(alpha))^2); % 角速度乘以比例系数
51
        tstart=te(end);             % 下一次求解的初始时间
52
        if abs(inity(2))<1e-6       % 如果此时角速度小于1e-6,则认为刚体已静止
53
            inity=[0,0];            % 将下一次的求解初始值设置为0
54
            nextIndex=find(abs(ag(index:end))>tan(alpha),1); %寻找地震波中下一个幅值超过tan(alpha)的点
55
            if ~isempty(nextIndex)
56
                index=index+nextIndex-1;
57
                sgn=sign(ag(index));
58
                tstart=interp1([ag(index-1),ag(index)],[tg(index-1),tg(index)],sgn*tan(alpha));  
59
            else                    % 如果找不到，则刚体不会再摆动，计算终止                
60
                ry(index,:)=ry_temp(end,:);
61
                break;
62
            end
63
        end
64
        continue;
65
    end
66
    if ie==2    % 如果由于摇摆角度超过alpha，则终止计算
67
        break;
68
    end
69
end
70
71
%% 微分方程
72
    function dy=wffc(t,y)
73
        dy=zeros(2,1); %用于存储计算结果：角度，角速度
74
        ag_t=interp1(tg,ag,t,'linear'); %根据线性插值，获取对应时间的地震加速度值
75
        %微分方程组
76
        %当刚体处于静止状态时，如果地震加速度小于tan(alpha)则不产生摆动
77
        if (abs(y(1))==0)&&(abs(y(2))==0)&&(abs(ag_t)<tan(alpha))
78
            dy(1)=0;
79
            dy(2)=0;
80
        else   %刚体摆动的微分方程
81
            sgn=sign(y(1));
82
            dy(1)=y(2);
83
            dy(2)=-p^2*(sin(alpha*sgn-y(1))+ag_t*cos(alpha*sgn-y(1)));
84
        end
85
    end
86
    
87
    % ode45的事件函数，当value=0时触发。
88
    function [value,isterminal,direction]=EventFcn(~,y)
89
        value(1)=y(1);      % 事件1：刚体摇摆方向改变，即此时刚体摇摆角度为0
90
        isterminal(1)=1;    % isterminal为1时，表示该事件的触发会导致求解终止
91
        direction(1)=0;     % direction为0，表示value上升到0，或者下降到0，都会触发事件
92
        
93
        value(2)=abs(y(1))-alpha;   %事件2：刚体摇摆角度超过alpha时，认为刚体倾翻
94
        isterminal(2)=1;
95
        direction(2)=1;             % direction为1，表示value上升到0时，才会触发事件
96
    end
97
end